
using JLD
using Plots
using Distributions

#-----------------------------------------------------------------------------
#DEFINE FUNCTIONS (cdfs and pdfs for Prob(m=L))
#-----------------------------------------------------------------------------
pr_fx(π,α;σ=0.011)  = Distributions.cdf(Normal(0,1),(α-π)/σ)
den_fx(π,α;σ=0.011) = Distributions.pdf(Normal(0,1 ),(α-π)/σ)
IMR_fx(π,α)         = den_fx(π,α) ./ pr_fx(π,α)
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
function Posterior_fx(α, Π, ζtil)
#-----------------------------------------------------------------------------
 zvL    = pr_fx(0, α).*ζtil ./ ( pr_fx(0, α).*ζtil + (1-ζtil).*pr_fx(Π , α) )
 zvNL   = (1-pr_fx(0, α)).*ζtil ./ ( (1-pr_fx(0, α)).*ζtil + (1-ζtil).*(1-pr_fx(Π , α)) )
 return zvL, zvNL
end
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
#Derivative of ζL wrt α
function Derivative_zetaL(Π,ζ,ζL, α,σ)
#-----------------------------------------------------------------------------
    num = IMR_fx(0.0,α) - IMR_fx(Π,α)
    den = pr_fx(0.0,α).*ζ + pr_fx(Π,α).*(1-ζ)
    der = (num./den) .* ( pr_fx(Π,α).*(1-ζ).*ζL ) ./σ
    return der
end
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
#Derivative of ζNL wrt α
function Derivative_zetaNL(Π,ζ,ζNL, α,σ)
#-----------------------------------------------------------------------------
    num = IMR_fx(-Π,-α) - IMR_fx(0.0,-α)
    den = pr_fx(0.,-α).*ζ + pr_fx(-Π,-α).*(1-ζ)
    der = (num./den) .* ( pr_fx(-Π,-α).*(1-ζ).*ζNL ) ./σ
return der
end
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
#Derivative of ζL, ζNL and ζL-ζNL wrt α
function Posterior_Derivative_fx(α, Π, ζtil)
#-----------------------------------------------------------------------------
#Compute posterior
zvL, zvNL  = Posterior_fx(α, Π, ζtil)
#Posterior derivative (for the normal CDF case)
der_zvL    =  Derivative_zetaL( Π, ζtil, zvL,  α, ce.σ_π)
der_zvNL   =  Derivative_zetaNL(Π, ζtil, zvNL, α, ce.σ_π)
der_zvL_NL = der_zvL-der_zvNL
return der_zvL, der_zvNL, der_zvL_NL
end
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
#Derivative of E(ζ(m)) wrt α: Surprise vs Frequency Channels
function Derivative_1st_fx(α,Π,ζL,ζNL,der_zvL, der_zvNL; σ_π=ce.σ_π) #wrt α
#-----------------------------------------------------------------------------
    Prob_Effect       =  (den_fx(Π,α) / σ_π) .* (ζL-ζNL)
    ProbL             =  pr_fx(Π,α)
    Surprise_Effect   =  ProbL .* (der_zvL - der_zvNL) .+ der_zvNL
    Surprise_Effect_L =  ProbL .* der_zvL
    Surprise_Effect_NL=  (1-ProbL) .* der_zvNL
    Total_Effect      =  Prob_Effect+Surprise_Effect
return Prob_Effect, Surprise_Effect, Surprise_Effect_L, Surprise_Effect_NL, Total_Effect
end
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
#LOAD BASELINE ECONOMY
#-----------------------------------------------------------------------------
α0 = -0.028
name_file = string("../01_Baseline_Model/model_data/CE_model_structure_alpha_",-α0,".jld");
ce = load(string(name_file), "ce");


#-----------------------------------------------------------------------------
#COMPUTE POSTERIORS & EXPECTED POSTERIORS FOR DIFFERENT {π,Π,α}
#-----------------------------------------------------------------------------
N1   = 100;
ζvec = [0.2 0.3 0.4 0.5 0.6 0.7];
πvec = linspace(-0.06,-0.0001,N1);
αvec = [-0.034; -0.030; -0.028; -0.026; -0.024; -0.022];
N2   = length(αvec);
N3   = length(ζvec);

Effect_Total_1st        = zeros(N3,N1,N2);
Effect_Message_1st      = zeros(N3,N1,N2);
Effect_Surprise_1st     = zeros(N3,N1,N2);
Effect_Surprise_1st_L   = zeros(N3,N1,N2);
Effect_Surprise_1st_NL  = zeros(N3,N1,N2);


for (iα,αv) in enumerate(αvec)
for (iπ,πv) in enumerate(πvec)
for (iΠ,Π) in enumerate(πvec)
for (iz,ζtil) in enumerate(ζvec)
		zvL, zvNL = Posterior_fx(αv, Π, ζtil)

        #Derivatives of the posteriors wrt α
        der_zvL, der_zvNL, der_PL_PNL  =  Posterior_Derivative_fx(αv, Π, ζtil);

        #1st derivative wrt α  (around π=Π)
        Prob_Effect, Surprise_Effect, Surprise_Effect_L, Surprise_Effect_NL, Total_Effect =
        Derivative_1st_fx(αv,Π, zvL,zvNL, der_zvL, der_zvNL);
        Effect_Surprise_1st[iz,iΠ,iα]    = Surprise_Effect
        Effect_Surprise_1st_L[iz,iΠ,iα]  = Surprise_Effect_L
        Effect_Surprise_1st_NL[iz,iΠ,iα] = Surprise_Effect_NL
        Effect_Message_1st[iz,iΠ,iα]     = Prob_Effect
        Effect_Total_1st[iz,iΠ,iα]       = Total_Effect
end
end
end
end

#-------------------------------------------------------------------
#-------------------------------------------------------------------
#PLOT
#-------------------------------------------------------------------
#-------------------------------------------------------------------
iz     = searchsortedfirst(ζvec[:],0.5)
iπ_min = searchsortedfirst(πvec,-0.050)
iπ_max = N1
#PLOT: 1st derivative
	iα = searchsortedfirst(αvec,α0)
	plot(πvec[iπ_min:iπ_max],  Effect_Surprise_1st[iz,iπ_min:iπ_max,iα]/100, label="Surprise Channel", linestyle=:dash, color="red")
	plot!(πvec[iπ_min:iπ_max], Effect_Message_1st[iz,iπ_min:iπ_max,iα]/100, label="Frequency Channel", linestyle=:dash, color="blue")
	plot!(πvec[iπ_min:iπ_max], Effect_Total_1st[iz,iπ_min:iπ_max,iα]/100, label="Total Effect", linestyle=:solid, color="black", linewidth=2.5)
	plot!([0], seriestype = :hline, color=:black, label="", legend=:topleft)
	ylabel!("Derivative")
	ylims!(-0.28,0.28)
    xlims!(-0.049,0.)
	xlabel!("Misreport Policy")
	xticks!([-0.05, -0.04, -0.03, -0.02, -0.01, 0.], ["-0.05", "-0.04", "-0.03", "-0.02", "-0.01", "0.0"])
	savefig("figures/Derivative_1stalpha.pdf")
#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------